uniform sampler3D noiseTex;
uniform sampler3D flowTex;

void euler();
void runge_kutta();

//to use, create a "standard" shader program:
//VertexWithTexture.vert +
//FlowVis3D.frag
//use 3D noise texture noise_tex_3D.raw (64x64x64)

void main()
{
	//euler();
	runge_kutta();
}

void euler()
{
	//min value for the tornado dataset
	//should be uniform var in general case
	float adj = -0.3;
	float stepSize = 1.0/48.0;

	vec4 result = vec4(0.0, 0.0, 0.0, 0.0);
	
	vec4 V = texture3D(flowTex, gl_TexCoord[0].stp); //velocity at [x,y,z]
	vec3 currentPosition = gl_TexCoord[0].stp;
	
	int num_steps = 50;

	for(int i = 0; i < num_steps; i++) //forward integration
	{		
		vec4 noiseTexResult = texture3D(noiseTex, currentPosition);
		if(noiseTexResult.a < 0.5)//noise texture is interpolated...
		{
			noiseTexResult.x = 0.0;
			noiseTexResult.y = 0.0;
			noiseTexResult.z = 0.0;
		}
		else
		{
			noiseTexResult.x = 1.0;
			noiseTexResult.y = 1.0;
			noiseTexResult.z = 1.0;
		}
		noiseTexResult.a = 1.0;
		result += noiseTexResult;
		
		V = texture3D(flowTex, currentPosition);
		V.x += adj;
		V.y += adj;
		V.z += adj;
		currentPosition.x += V.x * stepSize;
		currentPosition.y += V.y * stepSize;
		currentPosition.z += V.z * stepSize;
	}

	V = texture3D(flowTex, gl_TexCoord[0].stp);
	currentPosition = gl_TexCoord[0].stp;
	for(int i = 0; i < num_steps; i++) //backward integration
	{
		vec4 noiseTexResult = texture3D(noiseTex, currentPosition);
		if(noiseTexResult.a < 0.5)//noise texture is interpolated...
		{
			noiseTexResult.x = 0.0;
			noiseTexResult.y = 0.0;
			noiseTexResult.z = 0.0;
		}
		else
		{
			noiseTexResult.x = 1.0;
			noiseTexResult.y = 1.0;
			noiseTexResult.z = 1.0;
		}
		noiseTexResult.a = 1.0;
		result += noiseTexResult;
		
		V = texture3D(flowTex, currentPosition);
		V.x += adj;
		V.y += adj;
		V.z += adj;
		currentPosition.x += V.x * stepSize * (-1.0);
		currentPosition.y += V.y * stepSize * (-1.0);
		currentPosition.z += V.z * stepSize * (-1.0);
	}
	
	result.a *= 0.02;
	gl_FragColor = result / (float(num_steps)*2.0);
}


void runge_kutta()
{
	//min value for the tornado dataset
	//should be uniform var in general case
	float adj = -0.3;
	float stepSize = 1.0/48.0;

	vec4 result = vec4(0.0, 0.0, 0.0, 0.0);
	
	vec4 V = texture3D(flowTex, gl_TexCoord[0].stp); //velocity at [x,y,z]
	vec3 currentPosition = gl_TexCoord[0].stp;
	
	int num_steps = 40;

	for(int i = 0; i < num_steps; i++) //forward integration
	{		
		vec4 noiseTexResult = texture3D(noiseTex, currentPosition);
		if(noiseTexResult.a < 0.5)//noise texture is interpolated...
		{
			noiseTexResult.x = 0.0;
			noiseTexResult.y = 0.0;
			noiseTexResult.z = 0.0;
		}
		else
		{
			noiseTexResult.x = 1.0;
			noiseTexResult.y = 1.0;
			noiseTexResult.z = 1.0;
		}
		noiseTexResult.a = 1.0;
		result += noiseTexResult;
		
//		V = texture3D(flowTex, currentPosition);
//		V.x += adj; V.y += adj; V.z += adj;
		
//		currentPosition += V.xyz * stepSize;
		
		//
		vec3 k1 = V.xyz * stepSize;
		vec3 tempPosition = vec3(currentPosition.xyz);
		tempPosition += k1 / 2.0;
		V = texture3D(flowTex, tempPosition);
		V.x += adj; V.y += adj; V.z += adj;
		vec3 k2 = V.xyz * stepSize;
		
		tempPosition = currentPosition;
		tempPosition += k2 / 2.0;
		V = texture3D(flowTex, tempPosition);
		V.x += adj; V.y += adj; V.z += adj;
		vec3 k3 = V.xyz * stepSize;
		
		tempPosition = currentPosition;
		tempPosition += k3 / 2.0;
		V = texture3D(flowTex, tempPosition);
		V.x += adj; V.y += adj; V.z += adj;
		vec3 k4 = V.xyz * stepSize;
		
		currentPosition += k1/6.0 + k2/3.0 + k3/3.0 + k4/4.0;
	}
	
	result.a *= 0.02;
	gl_FragColor = result / (float(num_steps));
}

